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The topological response to external perturbations is an effective probe to characterize different 
topological phases of matter. Besides the Hall conductance, the Hall viscosity is another example of 
such a response that measures how electronic wave functions respond to changes in the underlying 
geometry. Topological (Chern) insulators are known to have a quantized Hall conductance. A 
natural question is how the Hall viscosity behaves for these materials. So far, most of studies on 
the Hall viscosity of Chern insulators have focused on the continuum limit. The presence of lattice 
breaks the continuous translational symmetry to a discrete group and this causes two complications: 
it introduces a new length scale associated with the lattice constant, and makes the momentum 
periodic. We develop two different methods of how to implement a lattice deformation: (1) a lattice 
distortion is encoded as a shift in the lattice momentum, and (2) a lattice deformation is treated 
microscopically in the gradient expansion of the hopping matrix elements. After establishing the 
method of deformation we can compute the Hall viscosity through a linear response (Kubo) formula. 
We examine these methods for three models: the Hofstadter model, the Chern insulator, and the 
surface of a 3D topological insulator. Our results in certain regimes of parameters, where the 
continuum limit is relevant, are in agreement with previous calculations. We also provide possible 
experimental signatures of the Hall viscosity by studying the phononic properties of a single crystal 
3D topological insulator. 


I. INTRODUCTION 

Topological phases are unusual phases of matter 
which cannot be categorized using Landau’s symmetry¬ 
breaking paradigm 1,2 . The recent discovery of topologi¬ 
cal insulators (TI) has opened new avenues for the study 
of symmetry-protected topological phases 3,4 , which re¬ 
quire some additional symmetries in order to remain in 
a robust topological phase. In all of these cases, the 
topological phase is protected by a bulk energy gap, 
e.g., in the case of the quantum Hall effects (QHE), the 
gap and non-trivial topology are due to an interplay of 
strong magnetic held and possibly electron-electron inter¬ 
actions. While topological phases have provided a fertile 
ground for theoretical study, there are also many exper¬ 
imental examples of various topological phases 5,6 . 

Since topological phases do not have a local order pa¬ 
rameter, one of the key goals for describing topological 
phases is to find a minimal set of quantities that fully 
characterizes each distinct phase. Such quantities pro¬ 
vide a framework to classify topological phases, and more 
importantly, may be related to the physical response co¬ 
efficients that can be measured experimentally, or used to 
distinguish topological phases in numerical simulations. 
A canonical example is the quantized Hall conductance 
for the QHE 5, 7 ~ 9 , while another example is the topologi¬ 
cal magneto-electric effect for 3D time-reversal invariant 
TIs 10,11 . 

Here, we are interested in another example of a re¬ 
sponse coefficient called the Hall viscosity 12 . This vis¬ 
cosity response is one part of the electronic viscoelastic 
stress response to an applied strain Uif 

\Tij ) AijkiUkl + Cij kl'U'kl j 

where T is the stress tensor, A and ( are the elas¬ 


ticity and viscosity tensors respectively, and the strain 
u ij = + djUi) is a symmetrized gradient of the 

displacement Ui . The Hall viscosity is one contribution 
to £ and is anti-symmetric in exchanging the first ij and 
second kl pairs of indices, and hence non-dissipative. For 
our purposes, we only consider isotropic systems (or sys¬ 
tems with at least C4 rotation symmetry), where the only 
non-zero components of the Hall-viscosity response are 
C1112 = C1222; which we denote by (A. Similar to the 
Hall conductance, <A can be derived from an adiabatic 
response calculation and is completely quantum in na¬ 
ture; it also requires time-reversal symmetry to be broken 
to be non-vanishing 12 . 

The Hall viscosity was originally studied in the integer 
and fractional QHE 12-41 and, for the integer effect, was 
shown to be <A = his 2 / 8 fI 2 b , where is = nh/eB is the fill¬ 
ing fraction, and Is = (fi/jeHl) 1 / 2 is the magnetic length. 
The viscoelastic response of the Chern insulator, a model 
for the integer QHE without Landau levels 42 , has also 
been investigated using a massive Dirac model 43 , and it 
has been shown that the Hall viscosity in the non-trivial 
Chern insulator phase is of the form C, = fi/ 8 tt£ 2 , where 
the emergent length scale depends on the bulk mass gap, 
to, and the Fermi velocity vf, as I = Hvf/Ziti. In both 
the integer QHE and Chern insulator one can write a 
continuum field theory from which the Hall viscosity is 
calculated as a response to deformations in the under¬ 
lying geometry of the system. In the former, the mag¬ 
netic length Ib introduces a minimal length/area scale, 
and (A is inherently finite; while in the latter, a regular¬ 
ization scheme was proposed to obtain the finite result 
above 43,44 . There have also been some recent studies of 
Hall viscosity response in 3D as well 45,46 . 

One open problem is the calculation of the Hall viscos¬ 
ity in lattice models away from the continuum limit. For 
example, in a discrete lattice tight-binding model it is not 
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clear, a priori , how to precisely model the geometric de¬ 
formations needed to generate a Hall viscosity response. 
In this article, we attack this problem by developing and 
comparing two different numerical methods to compute 
the Hall viscosity in tight-binding lattice models. These 
methods differ in the way geometric deformations are 
treated in the lattice model, and represent two natural 
lattice geometry modifications. Our first method is in¬ 
spired by the analogy between the viscoelastic and elec¬ 
tromagnetic responses, where we essentially introduce a 
frame field Peierls factor. In the second method, we use a 
more physical approach for how the strain can be realized 
by imagining mechanically moving the lattice degrees of 
freedom from their equilibrium positions and calculating 
the modifications to tight-binding overlap integrals. 

There have been some earlier pioneering works on Hall 
viscosity calculations in lattice models. First, the idea 
of coupling the electronic structure to phononic degrees 
of freedom, along the same lines as our second method, 
was originally studied by Barkeslili et al. 2A 1 where they 
proposed the “phonon Hall viscosity” as the adiabatic 
response of the electron state to acoustic phonons. They 
estimated the corrections to the phononic dispersion due 
to the viscosity, and found that it could be measurable in 
a number of materials, particularly ferromagnetic insula¬ 
tors. Other previous studies on lattice models 2 ' ,4 ' ,48 have 
calculated the Hall viscosity using entanglement proper¬ 
ties (momentum polarization), and have used the results 
to study the properties of the representative topologi¬ 
cal ground-state wave functions. Another recent study 49 
suggests that the Hall viscosity of a lattice Chern insu¬ 
lator may be related to a length scale associated with 
the Berry curvature at high symmetry points of the Bril- 
louin zone. This is motivated by the idea that the low 
energy effective theory of Chern insulator is dual to an 
effective theory in reciprocal space subject to an effective 
magnetic field (Berry curvature). Interestingly, the end 
result of this last work matches the continuum regular¬ 
ized result mentioned above. 

Our twofold focus here is different from the previous 
work. We first address how strain deformations can be 
modeled in tight-binding lattice models where continu¬ 
ous translational symmetry does not exist, and, second, 
to what extent the lattice results coincide with the field 
theory calculations. In fact, our aim is not only to pro¬ 
vide an alternative framework for numerical simulations, 
but also to relate the abstract notions in the field theory, 
such as geometric deformations, to the mechanical prop¬ 
erties of the crystal, such as phonons, as was also done 
in Ref. 24. This is a step toward providing direct signa¬ 
tures of the Hall viscosity in experiments on TI materials. 
In addition, we present a 3D generalization of these ap¬ 
proaches and use it to study the surface Hall viscosity of 
3D TIs with magnetic surface layers. 

To illustrate our methods we investigate three non¬ 
interacting models: (1) the Hofstadter model, i.e., the 
lattice version of the integer QHE, which primarily serves 
to benchmark our methods; (2) the Chern insulator lat¬ 


tice Dirac model, and (3) a 3D time-reversal invariant 
TI with magnetic layers on the surfaces. After the strain 
is implemented in the models, the Hall viscosity is com¬ 
puted using the Kubo formula 50 in terms of the corre¬ 
lation function of stress operators Tij (see Ref. 18 for a 
comprehensive discussion on this): 


C h = lim - -J2 
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where the system size is L x L and \o) denotes a single 
particle eigenstate of the Hamiltonian H\v) = E v \v). 


Our calculations confirm that for the Hofstadter model 
in a weak (compared to the lattice scale) magnetic field, 
the lattice calculations coincide with the continuum ex¬ 
pression. However, we start to see deviations from the 
continuum limit as the magnetic field is increased and 
lattice effects become more important. In the Chern in¬ 
sulator model, we find that in the vicinity of the criti¬ 
cal points where the topological phase transitions occur, 
the lattice results coincide with continuum limit predic¬ 
tions; however, away from the critical points, the value 
of the Hall viscosity is harder to determine and depends 
on our method of calculation. We also show that the mo¬ 
mentum polarization approach 2 ' -4 ‘ and our first method 
yield numerically identical results near the phase transi¬ 
tion points. For the 3D TI, we find that the Hall viscos¬ 
ity at the surface of the 3D TI can be fit as a quadratic 
polynomial in the surface mass gap. The coefficient of the 
quadratic term also has an interesting dependence on the 
bulk gap. Our results for both the 2D Chern insulator 
and the gapped surface of a lattice Dirac model share 
the common feature that the Hall viscosity is continu¬ 
ous as the phase boundaries are crossed, and eventually 
asymptotes to zero deep in the trivial phase. The Hall 
viscosity due to the non-trivial 3D TI electronic structure 
adds an anomalous term to the surface phonon dynam¬ 
ics, the effects of which can in turn be experimentally 
observed. Our estimates show that these effects are well 
within the precision of current experimental technologies. 
We briefly discuss possible experimental scenarios where 
these effects can be investigated. 


Our article is organized as follows: In Sec. II, we 
present the two methods to describe the effects of strain 
in the tight-binding models. The three subsequent 
Secs. Ill, IV, and V report the results for the three 
model Hamiltonians mentioned above, possible experi¬ 
mental signatures of the Hall viscosity are discussed at 
the end of Sec. V, and conclusions are outlined in Sec. VI. 
There is also an appendix, which contains some extra de¬ 
tails of derivations. 



II. TIGHT-BINDING MODELS IN THE 
PRESENCE OF STRAIN 
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Here, we introduce two different methods for model¬ 
ing strain on a lattice in order to couple the electrons 
to variations in the background geometry, see Fig. 1 for 
example. We will be considering discrete lattice mod¬ 
els represented by a generic multi-orbital tight-binding 
model 


H = ^ c\t r 


'C r ' i 


where cj. is a row vector of electron creation operators 
corresponding to the orbitals located in a unit cell at r, 
and t rr i is an overlap matrix between two sites at r and 




FIG. 1. Deformed lattice where the distortion tensor is given 
by u >22 = Gx. Note that bonds along x direction remain un¬ 
changed while the bonds in y direction are stretched uniformly 
as a function of x. 


A. Minimal Coupling as a Gauge Field 

For our first method we show that the strain in a tight- 
binding model can be modeled as a generalized Peierls 
substitution. Inspired by the analogy between the vis¬ 
coelastic formalism and conventional Maxwell electro¬ 
magnetism given in the field theoretical discussions of 
Refs. 28, 43, and 44, we can look for a generalized min¬ 
imal coupling in lattice models that must fulfill the fol¬ 
lowing requirements: 

i. The continuum limit must coincide with the held the¬ 
ory. 

ii. It must respect the 2ir periodicity of the lattice mo¬ 
menta; in other words, the Brillouin zone needs to be 
well-defined during the deformation process. 

iii. The distortion held Wij, to be dehned below, must 
be a continuous variable. 

The third point is crucial as we are planning to study the 
variations of the Hamiltonian with respect to infinitesi¬ 
mal strain deformations and hence derive the stress op¬ 
erators. 

Before we proceed, let us briefly review the aforemen¬ 
tioned analogy in a more held-theoretical language. In 
Cartan’s formulation of differential geometry in two spa¬ 
tial dimensions, a change in the metric is described by a 
set of local frame helds (vielbein) ej(x) where j = 1,2. 
The /i-th component of the frame held e.j (x) is denoted by 
e(( (x) (here /x = 1,2) in the local coordinate basis {d x , d y ) 
i.e. e -(x) = ej(x)9 x + e|(x)9 y . The metric tensor is de¬ 
hned in terms of the co-frame helds e J (x) = 

(dual to the frame helds, such that e l (e j) = <5() as 
g y „ = Sije^el, where 6 *• = is the Kronecker delta. 
Hence, in the absence of spatial curvature the torsion 
tensor is determined by T 3 = ( de 3 ) llv = e^d^e 3 , where 
the spin connection has been chosen to be zero using 


the gauge-freedom in the absence of curvature. As dis¬ 
cussed in Ref. 44, the matter held is coupled to the strain 
through its momentum 

Pj -A P^ = pj - p (2.1) 

where e j = Sj — wj, and Wj = djU M is the distortion 
(unsymmetrized strain) tensor. The shift in the mo¬ 
mentum is reminiscent of the standard minimal coupling 
(pj — eAj ) to the electromagnetic gauge held Af how¬ 
ever, here the “gauge held” Wj (overall four terms since 
j = 1,2 and p = 1 , 2 ) is multiplied by the momentum. 
So, the “charge” associated with this gauge held is indeed 
the momentum p y . 

Now consider the lattice Hamiltonian written in recip¬ 
rocal space H = ^k c kMk)ck- A uniform distortion, 
, (henceforth we only use lower indices for the dis- 
tortion/strain helds in the lattice models, we will also 
be slightly imprecise and refer to a distortion as a 
strain, and distinguish from the symmetrized version of 
the strain tensor by the symbol alone) is implemented as 
a shift in the lattice momenta: 


where a is the lattice constant. This choice of cou¬ 
pling clearly satishes all the requirements. In the contin¬ 
uum limit, a —>■ 0 , the above expression is simplified to 
ki —> ki — Wjikj which is consistent with Eq. (2.1). Note 
that there is a problem if one simply uses ki —>• ki — Wjikj 
as the minimal-coupling prescription for a lattice model 
because it does not satisfy the last two requirements si¬ 
multaneously, i.e., the second condition is not met for 
arbitrary Wjj, and if we force it to meet this condition, 
then Wjt has to be quantized and is no longer continuous 
(violation of the third condition). Another remark in the 
definition of Eq. (2.2) is that the conserved charge associ¬ 
ated with this gauge symmetry is the discrete momentum 
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operator 


ators 


Pi = Ya X^Px+a,- ^ C x+a 3 -Cx) 


(2.3) 


n-(9) _ 
ij 


SH (s) 


<5wi 


it;=0 


whose reciprocal-space representation is 


Pj = ~ sin ( fc i«) c k c k ■ 


We shall put the superscript ^ for all quantities used in 
connection to this generalized Peierls gauge-formalism. 
Using the transformation rule in Eq. (2.2), one can easily 

relate the stress operators in reciprocal space, = 

Ek c k^j 9 '(^) c k, the particle current operators J7i(k), 
and we find 


Now that we have specified the coupling to geometry 
let us describe the Hall viscosity response coefficient aris¬ 
ing from this method. In order to compute the Hall vis¬ 
cosity as a susceptibility in the linear response formalism, 
we need to determine the generalized force (stress) oper- 


= sm(kia) 
T^^k) = sin(/cia) 


8h{ k) 


Ski 
Sh{ k) 


w—0 


Sk 2 


w—0 


= sin(/eia) Ji(k), 
= sin(/cia) 


which describe the flow of momentum. From the defini¬ 
tion of the Hall viscosity, Eq. (1.1), one can write 


Ag) 


2 v (Q,k|T 1 ( i ?) (k)|^,k)(/3,k|r i ( 2 s) (k)|a,k) 
L 2 m ^ (E ak - E pk ) 2 

aGocc. 

/3(Eunocc. 


2 (a,k|sin(fcia)J'i(k)|/3,k)( / 9, k|sin(fcio) l / 2 (k)|a,k) 

T 2 m 4 " ^ ~ E ^ 2 

aGocc. 

/3(Eunocc. 


— -plm ^2 (sin 2 (A:ia) + sin 2 (^ 2 a)) 

k 

aGocc. 

/3Eunocc. 


(a,k|^i(k)|/3,k)^,k|J 2 (k)|a,k) 

(E a k — Ep k ) 2 


^ ^(sin 2 (/cia) + sin 2 (/c 2 a))-7 r (k) 


k 

k 


(2.4) 


where /i(k)|a, k) = U Q k|a, k), a refers to the band index, 
and H^^(k) is defined as the viscoelastic adiabatic cur¬ 
vature in the gauge coupling approach. In the third line, 
we have used the assumed C 4 symmetry of the system 
(invariance under k\ -f-> k 2 ) and write the expression in 
a symmetric fashion by including the contribution from 
[7 22 ,7 2 i]. This assumption is only necessary to simplify 
our discussion and we could relax the rotation symmetry 
without much extra difficulty. The usual adiabatic cur¬ 
vature associated with U(l) phase of the wave functions 
(Berry curvature) is denoted by .F(k): 

r ,u OT V- (a, k|j7 1 (k)|/3, k) ((3, k| t J 2 (k)|a, k) 

■ F(k,= - 2Im JL — — 

/3(Eunocc. 


If we use the the gauge invariant formula of the Berry 
curvature we have 

C h ] = YJ 2 5Z( sin2 fcl + sin2 fc 2)eijtr(T’ k 9 J T , k<9 J fPk) 

(2.5) 

where Uk = \ a ’ k)(a, k| is the projection operator 

onto the occupied states, and di = Using this re¬ 
lation will allow for a simple numerical computation of 

Ag) 
s h ■ 

Before we move on to the second method, it is useful to 
note that the existence of a conserved charge in this for¬ 
malism provides two alternative formulas for the Hall vis¬ 
cosity besides the linear response (Kubo) formula: (a) the 
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Streda formula and (b) momentum pumping/transport 
between the edges states using a Laughlin gauge argu¬ 
ment. Let us discuss these in a bit more detail. 

One can write the Streda formula for the Hall viscos¬ 
ity in terms of bulk quantities following analogy to the 
electromagnetic response. Let us motivate the idea by 
first reviewing the electromagnetic response. The Streda 
formula 51 for the Hall conductance is 


<Jh 


e 



T,B 


where e is the elementary charge, M z is the bulk mag¬ 
netization, [i is the chemical potential, T is temperature, 
and B denotes the magnetic held strength. The conju¬ 
gate quantities to the magnetization and the chemical 
potential are magnetic held and particle number density, 
respectively. By means of thermodynamic conjugacy re¬ 
lations we obtain 



Intuitively, this expression states that inserting a mag¬ 
netic hux binds charge to that flux. In a lattice model, 
one can look at the change in the particle density n(x) = 
E^occ \( X W(B))\ 2 as a function of the magnetic held 
(B) where the single particle states for a hxed B are de¬ 
noted by | v{B)). For small amplitudes of B the density 
varies linearly with B , and the Hall conductance can be 
read off from the slope. 

Using this as a guide, we expect that applying a vis¬ 
coelastic “magnetic held” would then add “momentum 
charge” into the system. The uniform viscoelastic “mag¬ 
netic held” G is determined by the viscoelastic “vector 
potentials” Wij. For instance, the viscoelastic vector po¬ 
tential/distortion tensor W 22 — Gx (which is in a Lan¬ 
dau gauge), corresponds to the lattice distortion shown 
in Fig. 1. For this choice of gauge, the system retains 
translational symmetry along the y-direction, and k 2 is a 
good quantum number. Therefore, the following Streda 
formula 28 can be proposed for the Hall viscosity 

c!?> = § (2.6) 

where the transverse momentum “charge” density is 
found by P 2 (x) = X^eocc. sin(/c 2 a)/a|(x|^(G))| 2 (x is in 
the bulk, far from edges) and |i/(G)) is the eigenstate of 
the Hamiltonian subject to the viscoelastic magnetic held 
of strength G. 

The second type of alternative viscosity calculation 
that can be used in the presence of a conserved charge 
relies on the momentum transport between opposite edge 
states when a viscoeleastic magnetic hux (dislocation) is 
threaded through the hole of a cylinder. To be explicit, 
let us consider edges realized in a cylindrical geometry 
(Fig. 2), where open boundary conditions are imposed 
along the x direction, and k 2 is still a good quantum num¬ 



FIG. 2. The generalized Laughlin experiment for the Hall 
viscosity. A viscoelastic flux is threaded through the cylinder, 
which one should think of as threading a dislocation through 
the cylinder. Electrons encircling the cylindrical hole will be 
translated by the Burgers’ vector of the dislocation flux. For 
the case illustrated here the Burgers’ vector is in the peri¬ 
odic y direction. If the system has a Hall viscosity there will 
be a momentum current flowing in the ^-direction carrying 
momentum pointing in the y-direction. 


ber in the periodic y direction. A generalized Laughlin 
experiment 8 can be done to measure the Hall viscosity. 
Threading a vicoelastic magnetic flux tpc along the axis 
of cylinder changes the transverse momentum argument 
k 2 of the Hamiltonian h(ki,k 2 ) into k 2 — (f>Gsin(k 2 a)/a. 
Thus, as 4>g is increased, the spectrum of the edge states 
is modified and there is a net momentum transfer from 
one end to another 44 . The Hall viscosity is given by 


(9) _ 1 d\P 2 
^ 2 d<t>G 


(2.7) 


in which A P 2 = P 2 ,r — P 2 ,l is the difference in the mo¬ 
mentum charge at the two edges of the cylinder. The 
momentum charge P 2 ,l(r) at the left (right) edge is cal¬ 
culated by summing the momentum density over a few 
unit cells (greater than the penetration depth of the edge 
modes) at the left (right) edge. The momentum density is 
given by P 2 (x) = J2^eocc. sin ( k 2 a )/ a \( x W(4'G)}\ 2 where 
\v(4>g)) is the eigenstate of the Hamiltonian subject to 
the viscoelastic magnetic flux of (\>g (Fig. 2). 

In Sec. IV, we will do an explicit calculation using all 
three of these sub-methods: Kubo, Streda, and Laughlin 
and show that they all match for the Chern insulator 
model. 


B. Electron Geometric Coupling Through Lattice 
Distortions 

For the second computational formalism, the strain is 
modeled based on microscopic deformations of a tight- 
binding model. This type of approach for the electron- 
phonon coupling in tight-binding models has been care¬ 
fully studied in graphene (see Ref. 52 and references 
therein) and there are interesting proposed effects based 
on the predictions of these calculations. For example, it 
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a 2 


1 


a i 



of the 3D TI when time-reversal is broken on the surface. 

Now let us introduce the method. A hopping matrix 
element t( r) in a tight-binding model is an overlap inte¬ 
gral between two orbitals spatially separated by r. We 
can define the strain field to be a function of the devia¬ 
tion dr from the equilibrium value ro- The linear order 
correction to the hopping matrix can be written as 


t{ r 0 + dr) « t( r 0 ) + 


(r 0 • dr) dt 


ro dr 


+ 0(dr 2 ), (2.8) 


FIG. 3. Illustration of the change in the overlap integrals as 
a result of strain. 


was shown that a non-uniform strain can lead to an effec¬ 
tive magnetic field with opposite signs at two valleys 53 , 
and this was subsequently experimentally observed 54 . 
This type of method was first applied to the Hall vis¬ 
cosity problem by Ref. 24. In that article they studied 
three models: (1) the lowest Landau band in the Hof- 
stadter model, although they ultimately only focus on 
the continuum limit; (2) a quantum spin Hall system 
in HgTe quantum wells 6 ’ 55-58 in the presence of a small 
time-reversal breaking magnetization, where they inves¬ 
tigate both the lattice and the continuum limit and dis¬ 
cover that there is a discontinuity in the Hall viscosity at 
the transition to the topological phase, and (3) a mean- 
field model of a p x + ip y superconductor 3,59 . A common 
assumption in this approach, which we will also use, is 
that in modeling the strain, the lattice distortion is con¬ 
sidered as an adiabatic process for electronic degrees of 
freedom. This implies that the phonon frequency must 
be much smaller than the electronic energy band gap. 
Moreover, we assume that the deformations are smooth 
on lattice scales, which implies the phonon frequency to 
be smaller than the Debye frequency. Note that the lat¬ 
ter assumption is not essential for lattice calculations and 
it is required only when one wants to take the continuum 
limit. 

In this subsection, we will derive a computational for¬ 
mula (Eq. (2.10)) for the Hall viscosity based on this 
model for the applied strain and we will apply it to two 
2D examples and one 3D example in subsequent sections: 
(1) the Hofstadter model where we study both numeri¬ 
cally and analytically the response of arbitrary integer 
filling fractions away from the continuum limit, and ( 2 ) 
the Chern insulator which is essentially half of the model 
considered for the quantum spin Hall effect in HgTe. In 
the latter case, however, we find there is no discontinuity 
at the transition to the topological phase, which we discus 
further below. We should emphasize that this observa¬ 
tion is also consistent with the field theory calculations. 
(3) We further obtain a 3D generalization of this formula 
and apply it to calculate the Hall viscosity at the surface 


which is due to the change in the bond length (isotropic 
contribution) and exists for all orbitals. In addition to 
this term, if more than one orbital (or local degree of 
freedom) is present in each unit cell, there can be a cor¬ 
rection to the hopping term between unlike orbitals due 
to an apparent rotation seen from neighboring sites. For 
instance, Fig. 3 shows a lattice model with two types of 
orbitals, s (blue) and p v (red), where the non-zero hop¬ 
ping terms along the 1 — 2 bond ai are only between alike 
orbitals. As a result of the deformation dr, the hopping 
matrix between s and p y orbitals becomes non-zero and 
proportional to the component of dr perpendicular to the 
unperturbed lattice vector ai, 

. . . n • (ai x dr) , 

ts,p y ( a i + dr) ss -—- t s , Py (a 2 ) 

l a il 

where n is the normal vector to the plane and f SjPy (a 2 ) 
is the hopping amplitude between s and p y orbitals in 
the vertical direction, a 2 . In general, dr is related to 
the strain tensor Uij (phonon field). For the example 
in Fig. 3, the hopping terms between sites 1 and 2 are 
modified as follows 


f s>s (ai + dr) w f S;S (a!) - u n t' s s 
tpy,Py( a 1 + ^ ~ ^Ps,P„( a l) — tpy,p y 
ts,p y (a i + dr) S3 U12 t StPy {a. 2 ) 

where = —adtg^/dr\ a is the derivative with respect 
to the lattice constant which is originally |ai| = |a 2 | = a. 
We remind the reader that we will assume C 4 symmetry 
for simplicity, and the overlap integrals could be modified 
in a more complicated manner if the lattice vectors are 
not orthogonal. 

The perturbed Hamiltonian can then be written as a 
function of and we define the strain operators by 


■fi p ) _ 
ij 


5&W 


5 Uij 


u —0 


on which the superscript ^ is placed for all quantities in 
this formalism. Following the linear response formula for 
the Hall viscosity, we start with 
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Ap) = 2 y- (t,a,k\Tii(k)\t,P,k)(t,P,k\Ti2(k)\t,a,k) 

^ L 2 m ^ (^ak(t) ^ ^k(i )) 2 

cnGocc. 

/3Eiinocc. 

2 y- (t,a,k|[a uil ,/it(k)]|t,/3,k)(t,^,k|[a ui 2 ,/i t (k)]|t,a,k) 

i 2 m ^ (£U(i) - ^k(i )) 2 

ccGocc. 

/3Gunocc. 

2 v (t,a,k|(a uil /i t (k) -/i t (k)a uil )|t,/3,k)(f, J 9,k|(a ttl 2 /i t (k) -/i t (k)5 Ul 2 )|t,a,k) 

L 2 m ^ (^ak(t) - ^/3k(0) 2 

a-Gocc. 

/3Gunocc. 

2 y^ (t,a,k|(L; ak - £ / 3 k )a uil )|t,/3,k)(i,/3,k|(£’ Qk - -E/ 3 k)< 9 Ul2 |t, a, k) 

£ 2 m 4 " (^k(t)-%W ) 2 

a-Gocc. 

/3Gunocc. 

= _ T2 Im Y (i ; a I k |3«n|t,,0, k )(t,/3, k |a„i 2 |^a, k ) 


k 

aGocc. 

ySGunocc. 

=^E B, ”(k) 

k 


where /i t (k)|t,a,k) = -Eak^lL cc, k) is the eigenstate of 
the Hamiltonian, h t ( k), with the set of hopping matrix 
elements collectively denoted by t. to labels the original 
(unperturbed) values of hopping amplitudes. The vis¬ 
coelastic adiabatic curvature B^ p \]s.) over the occupied 
states is also introduced. The above expression can be 
recast in a gauge invariant form 

C ( H = . (2.10) 

k 

Note that the dependence on comes from the depen¬ 
dence of the hopping matrix elements on the strain. We 
also urge the reader not to confuse the parameter depen¬ 
dence on the hopping amplitudes that we have denoted 
by t , and the time-coordinate, which does not enter any 
of the previous expressions. 


We have now completed the introduction of our meth¬ 
ods and in the following sections we will investigate sev¬ 
eral lattice models with continuum limits that have been 
a subject of great interest. The first model we look at is 
the Hofstadter model, which leads to the standard integer 
QHE Landau level problem in the continuum limit. The 
second and third models are described by massive Dirac 
Hamiltonians in the long-wavelength limit, i.e., the min¬ 
imal models for TIs. We compare the lattice results in 
each case with their continuum limit counterparts, and 
discuss the agreement between the two limits and why, 
or why not, we should have such an expectation. 


(2.9) 


III. EXAMPLE 1: THE HOFSTADTER MODEL 


Consider the Hamiltonian of a single band tight- 
binding model on a square lattice subject to a magnetic 
field 60 

h = -\Y telA ' j 4 \ Y ieiAiJ 4 c j ( 3 - 1 ) 

(i,J) «*,J» 

where { ) and (( )) refer to the nearest neighbor and next 
nearest neighbor sites respectively, and 

Aij = / A(x) • g?x 
J i 

with the choice of Landau gauge A = B(0,x). The mag¬ 
netic field is B = <^>/a 2 where the flux per plaquette is 
denoted by <f> = p/q ( p and q are coprime integers) in 
units of cj )o = h/e. We define a supercell of size q x 1 and 
derive the Hamiltonian in reciprocal space: 

H = Y C kM k ) c k 

k 

where c k = (c{ k ,..., c' q k ), and the first index labels the 
sublattice position within the magnetic supercell. The 
Hamiltonian matrix h(k) reads 













FIG. 4. The spectrum of the Hofstadter model for two values 
of the magnetic field (j>/<j>o = 1/20 (top) and 1/10 (bottom). 
Parameters are t = lOOt = 2, a = 1 and the system size is 
100 x 100. 


/ Ar(k) 3i(k) E q (k)e ikia \ 

Sr(k) A 2 (k) S 2 (k) 

S 2 (k) • 


s 9-t(k) 

\S g (k)e-^ a S,_!(k) A,(k) / 

in which A x (k) = —fcos(/c 2 a — Bax), and H x (k) = 
—1/2 — tcos(k 2 a — Bax — </>/2). The spectrum of this 
Hamiltonian is given by the Landau bands (see Fig. 4). 
We can then use Eqs. (2.5) and (2.10) to calculate the 
Hall viscosity. The results are plotted in Fig. 5. In the 
lowest Landau level, Fig. 5(top), there is a remarkable 
agreement with the continuum expression up to a field 
strength of 4>/<j>o = 3/50. After that, the effects of the 
lattice become prevalent and we find a deviation from the 
continuum expression. In Fig. 5(bottom), we compute (h 
for higher integer filling fractions and observe that the 
lattice calculations coincide with the continuum results 
up to v = 10. In both plots, the electron-phonon coupling 
method seems to stay closer to the continuum results over 
a wider range than the gauge coupling method. Since, a 
priori, we do not know what the correct value for the 
Hall viscosity is in a lattice system, we cannot identify 
which method, if either, is correct, only that they both 
reproduce the continuum limit, and both deviate when 
lattice effects become important. 




FIG. 5. Top: the Hall viscosity at v = 1 IQHE as a function 
of the magnetic flux per plaquette Bottom: the Hall 

viscosity vs the filling fraction v (integer fillings) where ^ = 
0.01 is kept fixed. C^h is in units of h/2na 2 where a is the 
lattice constant. 


Let us now show analytically that the continuum limits 
of both formulas lead to the same results. We will derive 
the continuum limit of the Hamiltonian in the absence 
of B, and then add B back in afterward. In the gauge 
coupling formalism, we expand the Hamiltonian around 
k = ( 0 , 0 ) as 


hgl 0 (k) + t)[(fci - wnkx - w 21 k 2 f 


+ {k 2 — w 22 k 2 — Wi 2 ki) 2 ) + ... 


1 


— ~ k-n^k 
“2m M « J 


where the mass is m = (t/2 +1 ) 1 / 2 , and the metric is 
given by 


n (g) = ( 1 ^ 2iuii -2wi2 \ 

y —2ui 2 i 1 - 2w 22 J 


(3.2) 


In the presence of the magnetic field, ki is simply replaced 
by ki—Ai (setting the charge e = 1). The stress operators 
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are found to be 


continuum limit in this formulation 


r (s) _ 

' n ~ m' 

T12 = „—( DiD 2 + D 2 Di), 


where Di = — — Ai is the electromagnetic covariant 

derivative. We introduce the ladder operators using the 
fact that [D 2 ,Di] = iB 


a = y ^g (£>2 +iDi), 

= V5 P2 “ W 1 ' )> 


and the strain operators become 


7 ^ ) = i ^(° 2 - 0ta ) 


2 to 

B 


(3.3) 

(3.4a) 

(3.4b) 


Let us compute the response for the n-th Landau level 


C H ] (n) 


2 T Y' ( n \Tn ) \n , )(n , \T 1 { 2 ) \n) 

2n£l ^ (E' n - E n ) 2 

B 2 y- (n|(a + a t ) 2 |n , )(n , |a 2 - a t 2 |n) 

(^ - E n y 

2n + 1 \ 1 

4 J 2tt£|' 


Here |n) is the n-th Landau eigenstate ada|n) = n|n), 
and the factor of 1 /L 2 was canceled by the degeneracy 
of the Landau level L 2 /2-k£ 2 b . Hence, for integer fillings 
v > 1, we can use the above results to derive 


d s) -5^d s) (n)- (3.5) 

n =0 & 

which is the same as the original findings of Avron et 
al. (and also Levay) 12-14 for the integer QHE viscosity 
response under modular deformations in the metric. 


Let us now proceed to the second method. The 
electron-phonon coupling formalism is modeled by the 
following changes in the hopping amplitudes 


d P io( k ) -i\ + *)( fc i + k l ) ~ + u 22 kj) 

t' 

— ^( M n + ^ 2 ) — 2t'ui 2 kik 2 


— ^ k n^k ■ 

~2m l9i i 3 


where the mass is defined as m = (t/2 +1) 1 /2 and the 
metric is 


Jp) 


= 1 - 


otiun + a 2 u 22 2a 2 u\ 2 

2a 2 ui 2 a\u 22 + 02*^11 


(3.6) 


where ai = m(t' + t'), and a 2 = mi'. Notice that the 
phonon field utj , as it appears in the metric, is multiplied 
by the electron-phonon coupling coefficients c^. We will 
see below that this will give a pre-factor in the final ex¬ 
pression for the Hall viscosity. 

Again, the effect of the magnetic field is recovered by 
substituting ki with ki — Ai. The stress operators can be 
found easily 

= + t')D 2 + t'D 2 ), 

rd = t'(DJJ 2 + D 2 D ± ), 

which can be rewritten in terms of the ladder operators 
(Eq. (3.3)) 

T 11 ^ + a d + ^ + at “))> 

Ti2^ = iBt'(a 2 — a^ 2 ). 


The Kubo formula for the n-th Landau level yields 

dV)=- 


*(p)^ _ 2 \ n ') (n 1 \Ti 2 ] \ri) 


2t:£ 2 b 


( K-En ) 2 


1 t't'B" ^ n{n — l)5 n ',n-2 


2 t t£ 2 b 2 


E 


( 2—) 2 

(n + l)(n + 2)5 rc/,n+2 


( 2^) 2 

k m' 


=t't’m 2 (2n + 1} 1 2 . 

4 2tt£ 2 


For filling fraction z/, the total Hall viscosity is found to 
be 


^x,x+ai ^ t ^11? 

^x,x+a 2 ^ ^ ^ ^ 22 ? 

^x,x+ai=ba2 ^ ^ ?(±«12 + ^(«11 +W 22 )) 

in which if = -off | r=a , t' = -\/2a|^| r=v ^ 0 , and we 
have used the identity M 12 = W 21 ■ Let us now look at the 


Ap) 


At't’ H 4 if if H 2 

(R^„ L „ 2n +1 = (§ + *) 2 8^ 2 / 

(3.7) 


where the usual Hall viscosity (2n + l)/87r^ is multi¬ 
plied by a non-universal dimensionless pre-factor which 
depends on the electron-phonon coupling constants if , 
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and i' and the band parameters. This does not contra¬ 
dict Eq. (3.5) since the phonon-induced changes in the 
underlying geometry of the electrons in Eq. (3.6) have 
extra coefficients depending on the electron-phonon cou¬ 
pling constants (as opposed to Eq. (3.2)), and hence the 
phonon Hall viscosity is proportional (and not equal) to 
the usual (gravitational) Hall viscosity. Note that above, 
when we showed the electron-phonon results in Fig. 5, 
we divided the data obtained from Eq. (2.10) by this 
pre-factor to compare only the magnetic field dependent 
factor. For the other models below we will not need to 
make this adjustment. 

We also note that the lattice effects on the Hall viscos¬ 
ity in the Hofstadter model has been studied recently in 
Ref. 61 using methods based on the momentum transport 
and momentum (entanglement) polarization. The Hall 
viscosity of Dirac electrons subject to a magnetic field 
has also been calculated 61,62 . In addition, a semi-classical 
derivation 29 for the Hall viscosity, which includes an anal¬ 
ysis of the Hofstadter model, has been worked out. It 
would be interesting, in future work, to carefully com¬ 
pare each approach. 


IV. EXAMPLE 2: CHERN INSULATOR MODEL 


zone. 

(b) Trivial phases: \m\ > 2r. The critical theory at 
the critical points between the trivial and topolog¬ 
ical phases are described by a single Dirac node at 
(0,0) or ( 7 r, 7 r) for m = 2 r and to = —2r respectively. 

For this model we calculate the Hall viscosity using both 
methods as to is varied. We show that the calculated 
viscosities have the following generic properties: 

i. C,h vanishes at the critical point to = 0 where the 
Cliern number changes from 1 to —1. 

ii. The Hall viscosity is continuous throughout the 
phase diagram. 

iii. The Hall viscosity is generically finite in the topo¬ 
logical phases and in the trivial phases, but becomes 
zero deep in the trivial phases (to —> ± 00 ). The 
rate at which (h goes to zero in the trivial phase is 
determined by the hopping amplitude t (which char¬ 
acterizes the gapless Dirac Fermi-velocity). 

All three properties are consistent with the field theory 
considerations in Refs. 43 and 44. 

Let us now discuss the details of the calculations. Us¬ 
ing the gauge coupling formalism, we write the deformed 
Hamiltonian in reciprocal space as 


For our next example we consider a (2+l)D lattice 
Dirac model for the Chern Insulator (Cl) on a square 
lattice 42,63 : 



E [ C x+a s (it s Vs 

X,S 


r<r 3 )c x + h.c. 


+ m Y 4ct 3 C x 

(4.1) 


where a s (s = 1,2) are the square lattice basis vectors, 
c x = ( C I xi c 2 x) is the electron creation operator with 
two orbital (or spin) indices, and (<Ti, er 2 , a 3 ) are Pauli 
matrices. The parameter r > 0 is the hopping amplitude 
between identical orbitals, and t\ and t 2 denote the hop¬ 
ping amplitudes between opposite orbitals (we choose the 
C 4 symmetric case ti = t 2 = t > 0). 

As the mass parameter to is varied, we obtain the fol¬ 
lowing phases: 


(a) Topological Chern insulator phases: — 2r < to. < 0 
with C = — 1 and 0 < m < 2r with (7 = 1. At the 
critical point between these phases (to = 0 ), there are 
two Dirac nodes at ( 0 , 7 r) and ( 7 r, 0) in the Brillouin 


=y^Ck ita 2 sin(A 2 - w 22 p 2 - w 32 pi) 


- ra 3 cos(A 2 - w 22 p 2 - w\ 2 pi) 


Ck 


Y c k * icr i s i n (^’i - W11P1 - w 21 p 2 ) 


- ra 3 cos(/ci - wnp 3 - w 21 p 2 ) 
+ TO E C^Ck 


Ck 


(4.2) 


where pi = sin A,; is the discrete momentum defined in 
Eq. (2.3). This expression immediately gives the stress- 

current operators T}f' > = Sk c k^ij 5 ^ (k)ck where 

T 11 (k) = {(Jit cos kia + a 3 r sin kia) sin Aqa, 

Ti 2 ^ (k) = (cr 2 t cos k 2 a + a 3 r sin k 2 a) sin Aqa. 

The lattice constant is denoted by a, and the total num¬ 
ber of sites is N = L 2 /a 2 . So, the Hall viscosity is given 
by (see Appendix A) 


d 9) = ^E s(s) ( k ) 

k 


1 

4iVa 2 


E 


f 2 (sin 2 k 3 a + sin 2 k 2 a){m cos Aqa cos k 2 a 


— r {cos Aqa + cos Aqa)) 
3/2 


[(to — r cos k 3 a — r cos k 2 a) 2 + t 2 (sin 2 A 3 a + sin 2 A 2 a)] 


(4.3) 
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Laughlin gauge argument calculations mentioned in Sec¬ 
tion II. In Fig. 7(a) we show the momentum density in 
a cylinder geometry as a function of the geometric de¬ 
formation for use in the generalized Laughlin experiment 
(c.f. Fig. 2). We see that the “momentum charge” with 
opposite signs is accumulated at the two edges of the 
cylinder as the viscoelastic magnetic flux 4>q is gradually 
dialed up. One can then read off the Hall viscosity using 
Eq. (2.7). We also compute by the Streda formula in 
Eq. (2.6). The outcomes of these calculations are plotted 
in Fig. 7(b), which shows that the Streda formula, as well 
as the Laughlin experiment, yield the same results as the 
Kubo formula. This remarkable agreement means that 
the concept of momentum charge in the gauge coupling 
formalism is still valid even in lattice models. 


Next, we derive the viscoelastic response using the 
electron-phonon coupling formalism. According to 
Eq. (2.8), the hopping matrix elements are modified such 
that 


t\(T\ -)• f( 1 - Ull)<Ti + tU21<J2, 
t 2 CT2 -t t( 1 - U22)P2 + tUi 2 ai, 
r x ->• r(l - Mn), 

r y -t r(l - U 22), 

where we have used the approximation dt(r)/dr ss 
— t/a 64 . Because of this simplification, and the form of 
the Dirac model coupled to strain, there will not be an 
extra pre-factor in the final expression of the Hall viscos¬ 
ity like what we saw for the Hofstadter model. From this 
we find the corresponding stress currents: 


FIG. 6. The Hall viscosity (in units of h/a 2 ) as a func¬ 
tion of mass for various values of hopping amplitude t in the 
range [0.25, 2]. The graphs are calculated from the Kubo for¬ 
mula in the thermodynamic limit N —> oo, (a) Eq. (4.3), (b) 
Eq. (4.4), and (c) momentum polarization method 27,47 . The 
arrow indicates the direction in which t increases with steps 
of 0.25. The colored regions show the topological phases and 
the Chern number ( C ) is indicated at the very top. The other 
hopping term is fixed r = 1 . 

We plot as a function of mass m for various val¬ 
ues of the hopping parameter t in Fig. 6 (a). Since in 
this formalism we have a conserved momentum charge, 
we can compare this result to the Streda formula and 


Tvi = \ E [4+ ai (^i - rcr 3)c x , 

X 

it \ r f + 

'12 =2 [ C x+a 2 ^l c x + c; +ai a 2 c x + h.c.J, 

X 

which can be written in reciprocal space as 

= cqfsinAqa — cos fcia, 

T(^(k) = t(cr 2 sin k±a + ay sin k 2 a). 

Plugging these expressions into the commutator of 
Eq. (1.1) leads to (see Appendix A for details of the 
derivation): 


= ^E B “( k ) 

k 


47Va 2 


E 


k 


t 2 (sin 2 A) 2 a (m — 2 r cos Ay a) + sin 2 k\a (m — 2 r cos k 2 Ci)) 
[(to — r cos k\a — r cos A^a ) 2 + f 2 (sin 2 k\a + sin 2 A^a )] 3 ^ 2 


Ap) 
S H 


(4.4) 
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4>G 




FIG. 7. (a) The illustration of the generalized Laughlin vis¬ 

coelastic gauge experiment (see Fig. 2). Each curve shows the 
momentum density profile (in units of h/a) along the cylin¬ 
der. The shift in the vertical direction is for presentation. 
Along the arrow from bottom to top the viscoelastic mag¬ 
netic flux cf>a is increased from (j>a = 0 to 4>g = 0.1 with steps 
of 0.01 . As a result, the momentum charge accumulated at 
the edges gradually increases. The bulk gap is fixed m = 1. 
(b) A comparison of the Hall viscosity (in units of h/a 2 ) in the 
gauge coupling formalism using the Kubo formula Eq. (4.3), 
the Streda formula, Eq. (2.6), and the spectral flow of the mo¬ 
mentum at the edges, Eq. (2.7). In both graphs, the system 
size is 100 x 100 and the hopping terms are r = t = 1. 


The results of the above expression are plotted in 
Fig. 6(b) as a function of m and t. As we see in this 
figure, there is no discontinuity around the critical point 
(|?7i| = 2r) from the topological to the trivial phases. 
This is in contrast with what was previously found in 
Ref. 24. Our calculation (details in Appendix A) shows 
that there is a factor of two for the ‘Vcos fc^a” terms in¬ 
side the parentheses in the numerator which is absent in 
the original calculation of Ref. 24. The difference in the 
results is attributed to this factor alone. If the factor 
is absent then the continuum limit of the Hall viscosity 
close to the phase transition (m —> 2 r) is different from 
Eq. (4.7) (as we will see) and the leading order would be 
linear in the mass gap, m — 2r. Inclusion of the factor 


of two, as our calculation indicates, yields results which 
match the regularized continuum model calculations of 
Refs. 43 and 44 that find that the lowest order of the 
mass gap entering the viscosity coefficient is quadratic. 

To gain a better understanding of what we have calcu¬ 
lated, we look at the Hamiltonian in the continuum limit. 
Let us start from Eq. (4.2). The electron operator can 
be written as 


Cx = ^ e iK,x^ (x) 

i 

where the slowly varying fields ^j(x), i = 0,1,2,3, are 
introduced around the four Dirac points at (0,0), (7r,0), 
(0,7r), and (7r,7r). The long-wavelength effective theory 
near each possible Dirac point can be written as 

#(c°nt.) = ^2 \°FOLi^(J V e^ v d^ + M,;cr 3 ]^i (4.5) 

where the Fermi velocity vf = ta , the masses are {Mi} = 
{m — 2r, m,m,m + 2r}, the sign coefficients are {a^i} = 
and {a iy2 } = and e? v is the 

frame field which can be written in terms of the distortion 
w 


^1,1/ — G^i^W 

to linear order in w. 

This is identical to the Hamiltonian considered in 
Ref. 43 using Pauli-Villars regularization. Therefore, the 
total Hall viscosity can be expanded as contributions 
from the neighborhood of each Dirac point 

c 

i 


where {CJ = { 0 ^ 10 ^ 2 } = {+, +} and 


I{M) 


k 2 M 


167T 2 


(vpk 2 


M 2 ) 3/2 


d 2 k 


M 2 M 2 + vik 


2 7,2 


87T 


Vpi/v 


|fc 2 


M 2 


M A M\M\ 

87 TVf 47 TVp 


(4.6) 


where a phenomenological cut-off A near each Dirac point 
has been introduced. It is worth mentioning that we have 
not set the sign coefficients C) arbitrarily as they appear 
in our calculations as a result of the long-wavelength ex¬ 
pansion around each Dirac point. Hence, 


j-(cont.) 
S H 



■ Amr) 


\m\ < 2r 
\m\ > 2r 


(4.7) 


in which the (arbitrary) A-dependent term has vanished. 

For comparison, if we take the continuum limit of 
the expression derived from electron-phonon coupling, 
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FIG. 9. Comparison of the Hall viscosity using the gauge 
formalism in Eq. (4.3) (solid) and the momentum polarization 
method” 1 , 47 (circles), t = 2.0 (blue) and t = 0.75 (red). 


FIG. 8 . The adiabatic curvature over the Brillouin zone in 
the vicinity of the critical point (m — 2 r = — 0 . 1 ). (a) and (c) 
are computed by Eq. (4.3) and (b) and (d) are computed by 
Eq. (4.4). In (a) and (b), t = 2 and in (c) and (d), t = 0.25. 
The energy scale is set by r = 1. The grid in reciprocal space 
is 200 x 200 . 


Eq. (4.4), this leads to the same results. For example, 
near the critical point, m = 2r + e, the difference in C,h 
between trivial and topological phases is 


AC£ ont ° 



(4.8) 


which is the same as the regularized field theory result 43 . 

A few remarks about Fig. 6 are in order. We have 
also included the results of momentum polarization 
method 27,47 in panel (c) for reference. The Hall viscosity 
is an odd function of the mass parameter ?n; this is consis¬ 
tent with the fact that the chirality is flipped as the mass 
changes sign for our model. However, the Hall viscosity 
does not vanish in the trivial phase as it does in the regu¬ 
larized field theory calculation. Interestingly, we do find 
that as the hopping coefficient t becomes smaller (i.e., 
Vp —> 0 ) (h goes to zero much more rapidly away from 
the critical point and into the trivial phase. Using these 
results we can interpret the discrepancy about the resid¬ 
ual viscosity in the trivial phase from several viewpoints. 
First, from a symmetry point of view, this result is not 
in contradiction with the requirements of generically hav¬ 
ing a non-zero Hall viscosity since time reversal symme¬ 
try is broken everywhere in the phase diagram. While 
the Chern number also requires time-reversal breaking, 
it does vanish in the trivial phase which is natural since 
the Chern number must be quantized, and is thus much 
more constrained. Moreover, crudely speaking, as we 
decrease t in the regime t < r, the effective bulk gap be¬ 
comes smaller and the corresponding correlation length 


becomes larger compared with the lattice constant. This 
essentially makes the lattice effects less important; con¬ 
sequently, we would expect (n to more closely match the 
continuum limit results where the Hall viscosity vanishes 
in the trivial phase. 

We can also understand the residual viscosity from an¬ 
other point of view by considering the viscoelastic adia¬ 
batic functions B^ 9 / p )(k) defined in Eqs. (4.3) and (4.4). 
For moderate regimes of t (when t ~ r) the quantities 
fi(g/p) (k) spread over a wider region of the Brillouin zone 
(Fig. 8 (a) and 8 (b)) while for smaller values of t, they 
become localized around the origin (Fig. 8 (c) and 8 (d)). 
Hence it is the latter case where we expect the contin¬ 
uum expansion to become more accurate, i.e., where this 
function is only sampling the contribution in the neigh¬ 
borhood of a Dirac point. 

Finally, from numerical point of view, the Hall viscos¬ 
ity is given in terms of the complex phase structure of 
the single-particle wave functions. If one of the three 
Pauli matrices was absent from the Hamiltonian, the 
wave function could be made completely real. In fact, 
deep in the trivial phases the 173 term is dominant over 
the entire Brillouin zone; hence all the wave functions are 
almost real in this regime and the Hall viscosity is negli¬ 
gible. In the vicinity of the transition point (|m| = 2r), 
but on the trivial side, the o\ and <72 terms can still be 
comparable with the <73 term within some regions of the 
Brillouin zone and this contributes to a non-zero Hall vis¬ 
cosity. As we tune down the hopping amplitude t, i.e., 
the amplitude of the o\ and ao terms, the region of the 
Brillouin zone over which these terms are comparable to 
the 03 term shrinks, and Qh asymptotes to zero faster. 
Hence because of the more complicated nature of the vis¬ 
cosity it is natural for it to be non-vanishing, even in the 
trivial phase. However, deep in the trivial phases, i.e., 
where we could consider the system in a trivial atomic 
limit, the viscosity indeed vanishes. 

It is interesting to note that besides the fact that the 
results of momentum polarization (Fig. 6 (c)) satisfy all 
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the aforementioned properties, it yields similar curves as 
the gauge formalism, a comparison is shown in Fig. 9. Re¬ 
markably, the difference between two methods are mini¬ 
mized as we approach the critical regions where the lat¬ 
tice effects are small. 


V. EXAMPLE 3: SURFACE STATES OF THE 3D 
TOPOLOGICAL INSULATOR 


For our final model we study the standard Wilson- 
Dirac Hamiltonian on a cubic lattice as a simple model 
of the 3D TI 10,65 

H =- ^2 4+a s {itsOLs - r/3)c x + h.c. 


s=l,2,3 

+ (to + 3r) ^2 4/3c x 


(5.1) 


where the Dirac matrices are given by 


Hamiltonian. For these phases it has been predicted that 
opening a gap in the surface states with a magnetic layer 
will induce a half-quantized quantum Hall effect 10,11 ’ 66 
and an accompanying surface Hall viscosity 43,45 . While 
the former prediction has been carefully studied in the 
context of the topological magneto-electric effect (also 
known as axion electrodynamics), the latter prediction 
has not been confirmed in a physical lattice model. 
Hence, our goal in this section is to explore the surface 
Hall viscosity response, and to also search for a route for 
experimental measurement by calculating the viscosity- 
modified surface phonon dispersion relations, as well as 
detecting non-local phononic responses (see Ref. 24 for 
similar phonon calculations for 2D systems). 

To calculate the viscosity we consider a slab geometry 
with a finite number of layers along the open-boundary 
^-direction, we keep the x and y directions translation 
invariant and periodic. We will induce a surface Hall vis¬ 
cosity by breaking TRS on the boundaries via TR break¬ 
ing Zeeman terms on the upper/lower surfaces 


a s = ri <g> a a 


0 <J S \ 

a s 0 J 


P = t 3 ® 1 = ( 

75 = Ti ® <7i = 


1 0 \ 

0 -1 J ’ 

'0 I\ 
V 1 0 )■ 


For our calculations we will take the Wilson parameter to 
always be fixed at r = 1 , and the system is taken to have 
cubic symmetry with ti = t 2 = t 3 = t. In this convention 
the <7 and r matrices act on the spin and orbital degrees 
of freedom respectively. 

Transforming to reciprocal space, the Bloch Hamilto¬ 
nian reads 


HpM — ^ O x ,z* 

x=(x,y) 

z=0,L z 

This term can be interpreted microscopically as layers 
of a ferromagnet deposited on the upper/lower surfaces 
that produce a magnetic field strength f l z . We take 
= — £Il z = H to have opposite signs, although this 
choice is not crucial, i.e., choosing the same sign yields 
identical results in our geometry though it might lead to 
complications in the presence of side surfaces. The fer¬ 
romagnetic term opens a Zeeman gap in the spectrum of 
topological surface states 


H u (i) = (Jiki =F a 2 k 2 ± (5.2) 


M k ) = 

3=1,2,3 


t s a s sin k s 


r/3 cos k s 


(to + 3 r)/3. 


There are two important symmetries present in this 
model:(i) time-reversal symmetry (TRS) with the oper¬ 
ator T = ia 2 IC such that a 2 h(k.)a 2 = h*{— k); (ii) inver¬ 
sion symmetry (IS), represented by I = t 3 V , such that 
T 3 h(\P)T 3 = h(— k). This model can exhibit a non-trivial 
3D TI phase protected by time-reversal symmetry. In 
fact, as the mass parameter m is varied, the Hamiltonian 
shows the following phases: 

(a) 0 < to and to < —6 r: trivial phase equivalent to the 
atomic limit. 

(b) — 2r < to < 0 and —6r < to < —4r: strong TI with 
a single Dirac cone on each boundary surface. 

(c) —4r < to < — 2r : weak TI with an even number of 
Dirac cones on each boundary surface. 

Here, we are only interested in the strong TIs in case 
(b) which have surface states described by a single Dirac 


where u(l) subscripts refer to upper(lower) surface states 
and are correlated with the signs. This Hamiltonian is 
valid in momentum space near the T-point up to a mo¬ 
mentum cut-off A, which depends on the bulk mass pa¬ 
rameter to, and the hopping parameter r. 

The gapped Dirac surface states on each surface have a 
Hall conductance of an = e 2 /2h which has been directly 
computed using several methods in Refs. 10, 11, 66 , and 
67. To calculate the viscosity, we can generalize one of the 
methods of Ref. 11, the “layer-resolved Chern number”, 
to compute the Hall viscosity layer by layer instead. The 
Hall conductance for a slab geometry can be written as 
an = Ce 2 /h, where C is the integer Chern number' given 

by 


C = 

k 

in which "Pk = Ylv l^kHiUkl is the projection operator 
onto the occupied states of the Hamiltonian in the slab 
geometry, and di = To find how different 2 layers 
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FIG. 10. Layer-resolved Hall viscosity using (a) the gauge 
coupling and (b) the electron-phonon coupling methods. Dif¬ 
ferent colors represent different values of the surface gap: 

= 0.0(green, triangles), 0.2(red, stars), and 0.4(blue, cir¬ 
cles). Other parameters are: m = — 1, r = t = 1, and 
(. L x ,L y ,L z ) = (100,100,20) 


contribute to C we can define a projection operator V z = 
\z)(z\ onto the c-th layer and compute 

2 t r 

C = ^2 E T r['PktiAW'P*(dj'Pu)\ ■ 

The same idea can be applied to the Hall viscosity. 
Starting from Eq. (2.5), the gauge coupling formula, one 
can calculate the contribution of layer z to the Hall vis¬ 
cosity as 

C h\z) = 2^-a Y ( sin2 fc i + sin2 fc 2) 

k 

x Tr[Vi l e ij (diP k )V z {d j 'P k )]. (5.3) 

Similarly, one can write down a layer-resolved expression 
for the Hall viscosity based on the electron-phonon cou¬ 
pling formula in Eq. (2.10), 

<h\*) = ^Y ^['P^(9 Uli V^)V z (d U2j V k )]. (5.4) 

k 

Figure 10 shows the layer-resolved Hall viscosity along 
the open boundary ^-direction using both methods. 
Next, the total Hall viscosity is calculated by summing 




FIG. 11. The Hall viscosity of the surface states of 3D TI. 
(a) The gauge coupling method Eq. (5.3). (b) The electron- 
phonon coupling method Eq. (5.4). The solid lines are fits 
using Eq. (5.5). The legend shows the bulk mass m values. 
The hopping parameters are t = r = 1. The system size is 
(/.,, /-„, l-r.) = ( 100 , 100 , 20 ). 


over the contributions from half of the lattice. 

The results of both methods for various values of f 1 are 
shown in Fig. 11. As we see, regardless of the bulk mass 
parameter to, the Hall viscosity is zero when the surface 
states are gapless Q = 0. Inspired by the continuum limit 
expression for the single Dirac cone in Eq (4.6), we can 
fit the data with the following ansatz 

<,5) 

where the coefficients (£i, £ 2 ) depend on the bulk TI gap. 
The sign change in curvature (sign of the quadratic term) 
as we pass from H > 0 to H < 0 is confirmed by our cal¬ 
culations. As shown in Fig. 12, both coefficients (^ 1 ,^ 2 ) 
decrease as the bulk mass is swept towards the trivial 
phase. From Eq. (4.6), and the fact that the cutoff is 
proportional to the bulk mass, the linear term should 
vanish as the bulk gap closes (to —»• 0 ); while as we see in 
Fig. 12(a), £1 stays non-zero beyond to = 0 and within 
the bulk trivial phase. This effect is similar to the one 
observed in 2D case where the Hall viscosity gradually 
vanishes in the trivial phase. We check that decreasing 
t indeed makes £1 transition to zero faster. In addition, 
the regularized 2D continuum expression in Eq. (4.6) im- 
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m 


FIG. 12. The coefficients of the fitting expression, Eq. (5.5), 
for the Hall viscosity in Fig. 11. Notice the coefficient of the 
quadratic term (£ 2 ) is much more method-independent than 
the linear term (£ 1 ). Parameter details at Fig. 11. 


plies that the quadratic term should not depend on the 
cutoff; therefore, we might expect that £2 = 1 and should 
remain constant in the topological phase. However, our 
lattice calculations in Fig. 12(b) show that £2 is close to 
1 only at the deepest point in the topological phase, and 
decays to zero as the bulk critical point is approached. 
We find that the coefficients of the quadratic term from 
both methods show similar magnitudes and trends, while 
the behavior of the nominally cut-off dependent coeffi¬ 
cient £1 (m) seems quite method dependent as we might 
expect from the continuum calculations. 

Now that we have shown that a Hall viscosity can ap¬ 
pear at the surface of a 3D TI, we will use the electron- 
phonon coupling formalism to discuss that the Hall vis¬ 
cosity can add new plrononic properties to the TI, which 
could be experimentally measurable. To this end, let us 
start with the effective theory for the phonon field. The 
phonon dynamics in the harmonic limit are governed by 
the following Lagrangian density 

£ph = 7; {^PUi CijklU-ijUkl^j , 

where p is the mass density and, for simplicity, we con¬ 
sider an intrinsically isotropic system where the elasticity 
tensor takes the form c ijki = X6 i:j S k i + p(S ik 6ji + 5u5 jk ). 
The parameters A and p, are called Lame constants 68 . 


According to the linear response theory, Eq. (1.1), the 
extra term due to the Hall viscosity is 

"^visc — 0/Gy U n iU n j (5.6) 

for a 2D Chern Insulator and, 

-^visc U n i3jU nk (h*7) 

for a 3D TI, which is an analog of the magneto-electric 
term E ■ B. 

We shall discuss the 2D case (which was originally ex¬ 
plained in Ref. 24) as a warm up example. The equation 
of motion is found to be 

Ui =v 2 V 2 u z + (v 2 - v 2 )diV ■ u + eijCn^Uj/p (5.8) 

where the velocities of the longitudinal (LA) and trans¬ 
verse (TA) acoustic phonons are defined via v( = 
a/(A + 2 p)/p and v s = \f~pjp, respectively. The Hall 
viscosity term couples the LA and TA modes and modi¬ 
fies the phonon dispersion relation as 

w « v a k + c a f —^ (— -—) k, (5.9) 

\ujh J v a 

for different modes a = £ or s, where the characteris¬ 
tic frequency determined by the Hall viscosity is u>h = 
p{v | — v 2 )/(h, and Q( s ) = ±1. The other consequence of 
the Hall viscosity term is that if we drive mode a with 
amplitude A a , the other mode a will acquire a non-zero 
relative amplitude of 

A a . u> 

A a U)H 

Notice that the Hall viscosity generated phonon wave has 
a 7r/2 phase shift which may make measurements easier. 

Now let us consider a 3D TI. For our 3D system let 
us assume a slab geometry with finite length along the 
^-direction. The equation of motion reads 

Ui =v 2 V 2 Ui + (vj - v 2 )di\7 ■ u + e ijk (dj(H)V 2 u k /p. 

(5.10) 

Note that unlike Eq. (5.8), the Hall viscosity term here, 
VC//(z) = C,h S(z)z, is a boundary term that modifies the 
surface phonon dispersion, as well as the boundary con¬ 
ditions. The correction to the dispersion relation (now 
only at the surface) has the same form as Eq. (5.9) up 
to replacing the 2D density p with the 3D density mul¬ 
tiplied by the penetration depth of the surface states 
(roughly few layers in the ^-direction, e.g. a few nm 

in Bi 2 Se 3 ), p —>■ where £ ~ hvp/E ( 'g Ulk ^ > is the sur¬ 
face state penetration depth. In addition, we discuss two 
unique signatures of the Hall viscosity in a 3D geometry: 
the analog of the Faraday rotation for normally incident 
TA phonons, and the emergence of a TA phonon wave 
traveling away from the surface due to the excitation of 
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surface LA phonons. 

First let us consider the phonon Faraday rotation. The 
general form of TA phonons propagating along the open 
boundary ^-direction is given by the displacement field 
u(z) = Ui(z)x + U2{z)y. Suppose TA phonons in the 
^-direction are excited, ui(z) = The 

boundary condition for the other TA mode is 


d z u 2 


2 = 0 + 


-i—^-rUJ 3 Ul(z = 0). 

pv t 


Note that in absence of the Hall viscosity, this equation 
reduces to the boundary condition for the free surface: 
d z u 2 = 0 and U 2 remains zero. However, for C,h ^ 0, this 
mode is generated, 112 = J 4 2 ed w /' u »)«- iw ^ with a relative 
amplitude of 


^2 C h 2 

- = — o w . 

Ai pv/ 

Hence the polarization of the phonon is rotated by a small 
angle 6 given by tan# = A 2 /A 1 in the ccy-plane. The 
rotation of oscillation axis is an analog of the Faraday 
rotation for the TA phonons. 

Next, let us now illustrate how boundary effects can 
lead to a non-trivial phononic effect (which is unique to 
phonons and does not exist for the photonic response) 
for a 3D TI. Consider the following general form for the 
phonon field u(x, z) = ui (x, z)x +112 (x, z)y . Suppose the 
surface phonons are driven in the LA mode as u\[x, z = 
0) = Thg boundary condition for the TA 

mode can be derived as 


d z u 2 


2=0+ 


■ C H 

PV 2 s V l 


uj 3 ui(z = 0 ). 


This equation implies that there will be a phonon wave 
traveling away from the surface, given by the ansatz 
u 2 (x,z) = A 2 e i ^/ Ve '> x+ik3Z ~ iut , where k 3 = - 

v^/v S V(. can be found from the wave equation in the 
bulk. The amplitude of such a mode can be computed 
from the boundary condition 


A 2 

Ai 


p(vj 


C H _ 

v%y/ 2 v s ve 


to 2 . 


This phenomenon could be considered as a direct signa¬ 
ture of the Hall viscosity in a 3D TI: the TA phonons 
can be detected on the bottom surface as a result of ex¬ 
citing the LA phonons on the top surface (illustrated in 
Fig. 13). 

To be a bit more quantitative we estimate the correc¬ 
tion due to the Hall viscosity as follows: for 2D systems, 
a typical solid has a mass density of p ~ 10 -7 g/cm 2 , 
phonon velocities of vi ~ 2v s ~ 10 5 cm/s, and C,h ~ 
I0~ 1 h/a 2 as we have found in the previous section. Con¬ 
sequently, the characteristic frequency is approximately 
ojh ~ 10 15 s _1 . The important assumption in our cal¬ 
culation is to treat the phonon in the adiabatic limit; 



FIG. 13. Traveling TA phonon wave into the bulk TI as 
a result of driving the surface LA phonons. The green and 
yellow stripes are to represent compression and decompression 
regions respectively. Arrows show the propagation directions. 


i.e., w -C ujD-,Eg, where a generic value for the Debye 
frequency is +>+> ~ 10 14 s _1 , and the electronic bulk gap is 
Eg ~ 10 14 s _1 . Hence, one can reach rather considerable 
ratios up to ui/uih ~ 10~ 2 at frequencies around 1 THz. 
For 3D systems, p ~ 10 g/cm 3 and (h ~ 10 _2 fi/a 2 as 
we have calculated in this section; therefore, the ratio 
A 2 /A 1 can be made as large as 10 -2 . Similar orders of 
magnitude, 10 -2 , can be accessed for the corrections to 
the surface phonon dispersion of 3D TIs. 

There are already various experimental techniques that 
might be used to investigate the unusual phononic prop¬ 
erties of TIs due to the Hall viscosity. In terms of ma¬ 
terials to be studied, Bi 2 Se 3 as an archetypal 3D strong 
TI could be a promising choice, especially since recent 
experiments have shown that the electron-phonon cou¬ 
pling at its surface is strong 69 . To study 2D effects (wave 
mixing), the surface acoustic (Rayleigh) waves method 
can be used to excite and measure different phononic 
modes. This method has been applied successfully to 
measure the electron-phonon coupling of a 2D electron 
gas (GaAs heterostructures) in integer QHE 70 and frac¬ 
tional QHE' 1 regimes. Recently, a similar setup was ex¬ 
amined for the study of phonons coupled to Dirac elec¬ 
trons in graphene' 2 . 

To explore the 3D properties (analog of Faraday ro¬ 
tation, and emergence of a normal propagating mode), 
there are pump-probe measurements where LA/TA 
phononic modes are excited as a result of the relaxation 
of hot carriers photoexcited with ultrashort (~ 100 fs) 
laser pulses. For the slab geometry, normally incident 
phonons are directly measured using different techniques 
such as time-resolved x-ray diffraction' 5,14 and transient 
reflectivity (TR) traces ,5 ~ 77 . The first method can be 
used to capture both LA and TA phonons while the sec¬ 
ond method is particularly useful for measuring the LA 
phonons. Recent experiments on phonon dynamics in 
Bi 2 Se 3 " can in principle be readily adapted to study 
the Hall viscosity effects by adding ferromagnetic sur¬ 
face layers or applying an external magnetic field. As 
a final remark, we note that crystal imperfections may 
also lead to mixing between LA and TA phonons, how¬ 
ever, they can be distinguished from the Hall viscosity 
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related effects because they do not depend on the sign of 
the time-reversal symmetry breaking terms such as the 
direction of magnetization in ferromagnetic layers or ex¬ 
ternal magnetic field. 

In summary, from these calculations we can expect a 
non-zero surface Hall viscosity, but its exact value and 
dependence on the bulk mass m is potentially quite com¬ 
plicated, and model dependent. We see some expected 
trends in the fitting coefficients that qualitatively match 
a continuum model for the surface states, and the proper 
symmetry constraints are met, but the results seem de¬ 
pendent on the choice of method, and on the high-energy 
regulator of the surface theory, i.e. the bulk mass to. 
The presence of a non-zero Hall viscosity, regardless of 
the complicated dependence on details of the band struc¬ 
ture, leads to various measurable effects in the phononic 
response of TIs: it modifies the phonon dispersion rela¬ 
tion (bulk phonons in case of a 2D Chern insulator, and 
surface phonons in case of 3D) and it mixes orthogonal 
phonon modes which would be otherwise independent. 
There are other interesting questions one may ask about 
the emergent Hall viscosity on the boundaries of 3D TIs. 
Here, we computed the surface contributions due to the 
Zeeman term imposing the periodic boundary conditions 
on the boundaries of the top and bottom surfaces in a slab 
geometry; one can remove this boundary condition to see 
how gapless chiral domain wall states on side surfaces 
contribute to the Hall viscosity when there is a magneti¬ 
zation domain wall. Another possibility is that one can 
include orbital magnetic held terms, as a result of which 
Landau orbits will form in presence of a strong magnetic 
held. We speculate that the massive Dirac surface states 
would have an anomalous Hall viscosity due to the non¬ 
zero Dirac mass, in addition to the usual Hall viscosity 
coming from the Landau orbits. As a slightly different di¬ 
rection, one can also study Hall viscosity purely due to 3D 
bulk when time-reversal is explicitly broken in the bulk. 
We expect that in this case viscoelastic response would 
be non-zero without any magnetic layers on boundary 
surfaces. 


VI. DISCUSSION 

We have presented two different schemes to imple¬ 
ment lattice deformations in order to calculate the Hall 
viscosity in tight-binding models. We have validated 
these methods showing that, numerically and analyti¬ 
cally, they both converge to the same results in the weak 
held limit of integer quantum Hall on the lattice (Hof- 
stadter model). Next, we have applied them to the 2D 
Chern insulator model and have demonstrated that, al¬ 
though different methods give rise to different values for 
the Hall viscosity, they satisfy certain common properties 
due to symmetry considerations; i.e. the Hall viscosity 
hips sign as the Chern number changes sign, vanishes at 
the particle-hole symmetric point to = 0 where the Chern 
number changes sign, is continuous across the topologi¬ 


cal phase transitions, and asymptotes to zero in trivial 
phases. A comparison with the momentum polarization 
method shows that the gauge formalism matches well 
with the momentum polarization results near the phase 
transition points where the lattice effects are minimized. 
We illustrated that confining the Berry curvature to the 
origin of Brillouin zone k = 0 will force the lattice Hall 
viscosity to look more like the continuum Hall viscosity 
calculated for massive Dirac fermions in 2D. In addition, 
starting from the lattice tight-binding model, the contin¬ 
uum limit is derived and shown to yield similar results 
to previous continuum limit studies 44 . We have gener¬ 
alized our methods to 3D models and, in particular, we 
have studied the surface Hall viscosity of 3D TI where 
the surface Dirac states are gapped as a result of a Zee- 
man splitting term. Inspired by the expression for the 
Hall viscosity of single massive 2D Dirac fermion in the 
continuum, we fit our data and find that the nominally 
cut-off independent term is much more method indepen¬ 
dent as well. The overall feature of our calculations in 
two or three dimensions is that the Hall viscosity, unlike 
the Hall conductance, is continuous as one crosses the 
critical points between topological and trivial phases (or 
from topological to topological phase). 

It is worth noting that the difference between the re¬ 
sults of our two methods away from critical points are 
most likely due to lattice effects and not related to the 
possible differences between a conventional Hall viscosity 
and a torsional Hall viscosity as discussed in Refs. 33 and 
34. The conventional Hall viscosity refers to those calcu¬ 
lations in systems which usually have Galilean invariance 
(including e.g., rotational invariant QHE states and chi¬ 
ral superfluids), and the torsional Hall viscosity refers 
to those calculations in the Dirac-type models with cou¬ 
pling between the momentum and spin/orbital matrices. 
The latter requires a coupling to geometry through the 
frame-fields, while most aspects of the former are satis¬ 
factorily treated using just the metric tensor. Both types 
of calculations share many similar features, but both ap¬ 
proaches have yet to be unified into a full understanding. 
In this respect, our calculations (since they reproduce 
the correct continuum-limit values of both of the poten¬ 
tial types of Hall viscosity) suggest that there should be 
no distinction between the torsional and the conventional 
Hall viscosities, at least as far as the stress-stress response 
is concerned. This idea is supported by the observation 
that the continuum limit of our two methods (along with 
the third method based on the momentum polarization) 
in the Hofstadter or the Chern insulator models lead to 
the identical correct results in both cases. Another piece 
of evidence is based on our calculations for the surface 
of a 3D TI, where the Hall viscosity computed by the 
two methods has an almost method independent piece 
(quadratic in surface gap). 

We note that the underlying idea for our derivation 
of the Hall viscosity is quite general and does not in¬ 
volve any assumption about the type of lattice, crystal 
symmetry group and the shape of sample. Therefore, 


19 


our techniques can be adapted to wide variety of lattice 
configurations and sample geometries. Remarkably, in 
the presence of translational symmetry, one may write 
the Hall viscosity in terms of an integration over lattice 
momenta which can be evaluated exactly in the thermo¬ 
dynamic limit by means of well-known numerical integra¬ 
tion techniques. 

Using the electron-phonon coupling (second) approach, 
we have shown that the Hall viscosity would appear as an 
anomalous term in the effective action for phonons. This 
term couples the longitudinal and transverse phononic 
modes and modifies the dispersion relation. We have ex¬ 
emplified various possible experimental scenarios to ob¬ 
serve direct signatures of the Hall viscosity by investigat¬ 
ing the phononic properties of a media. 

In principle, our methods can be applied to all other 
interesting gapped or gapless lattice models for interact¬ 
ing or non-interacting topological phases. Both formulas, 
Eqs. (2.5) and (2.10), for the Hall viscosity can be gen¬ 
eralized to the many-body problems 

C H ] = -2 Im (<9 Ul2 ’I'|<9„ 11 ’F) 

where I’l') is the many-body wave function. 

Last but not least, provided that the Hall viscosity can 
be used as another probe along with other topological 
responses to distinguish topological phases; one can also 
study the effect of interactions (and maybe disorder) on 
the viscoelastic behavior using the methods introduced 
here. 


Here for simplicity, we set (the lattice constant) a = 1. 
For each k, there are two eigenstates 


H\±,k)=±E k \±,k), 


where 
I- k) = 

and 


w k 

-u k e 


i6 k I ! 


■ k) = 


Vk 

u k e l6k 


Ek=(h 2 1 (k) + h 2 2 (k) + h 2 3 (k)) 1/2 , 
IE k - h 3 { k) 


Uk 


Vk = 


2 £ k 


I E k + h 3 (k) 


2E k 


tan 9 k = 


Mk) 

hi(k)' 


(A3) 


The matrix elements of the Pauli matrices in the eigen¬ 
state basis are 

-+ / ii i, i\ ?’£k sin d k — h 3 cos 9 k 

oi = (-,k|<7i|+,k) =- - -, (A4a) 

E k 

-+ / ii i, i\ -LE k cos0 k - /i 3 sin0 k 

= <-,k|cr 2 |+,k) = ---, (A4b) 


Ek 


< J 3 + = (-,k|<T 3 |+,k) = 


Ek 


(A4c) 


VII. ACKNOWLEDGMENTS 


1. Minimal coupling method 


HS would like to thank Maissam Barkeshli, Pouyan 
Ghaemi and Thomas Tuegel for useful discussions. TLH 
is supported by NSF (USA), DMR-1351895-CAR. SR has 
been supported by the grant NSF (USA), DMR-1455296 
and Alfred P. Sloan foundation. 


We are to calculate 

(-,k|7n s) |+,k) = (d 1 h 1 a^ + + d 3 h 3 a 3 + ) sin ki, 
(—, k|7i2^|+> k) = {d 2 h 2 cF 2 + + d 2 h 3 a 3 + ) sinfci. 


Appendix A: Derivation of the Hall viscosity 
formula for the Chern Insulator 


In order to compute the Hall viscosity using Eq. (1.1), 
we need to find the eigenstates of the Hamiltonian, 
Eq. (4.1). The Hamiltonian in the reciprocal space reads 

H = cj c /i s (k)<T ; ,.c k (Al) 

s, k 


where 


Nk =(-, k|rA 9) |+, k)(+, k\T^ I-, k) 
d\hid 2 h 2 

- ^2 - y l ih 2 + lh 3 Ek) 

d\h 3 d 2 h 3 2 , ,2\ 

+ -^2- ( h l+ h 2) 

V'k 

dihid 2 h 3 

H-^-( iEkh 2 — h 3 h i) 


d\h 3 d 2 h 2 


(iEkhi — h 3 h 2 ) sin 2 fci 


hi{ k) = t 3 sin ki, 
h 2 (k) = t 2 sin k 2 , 
h 3 ( k) = m — r cos k± — r cos k 2 . 


and the Hall viscosity becomes 


(A2) 


d 9) = ^E 6(9) ( k ) 

k 


(A5) 


(A6) 
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where the adiabatic curvature is 

B<S,(k >=“ 2 I "'{( 2 t?} 

E h a dihbd2h c . 2 , 
£ahc - wf * - sm «i 


abc 


2 K 


tit 2 sin 2 ki (to cos k\ cos fc 2 — r cos k\ — r cos fc 2 ) 


2 K 


and then we can make it isotropic 
B (a)( k ) - t i^(sin 2 fc^+sin 2 fc 2 ) ;; 


(to cos fci cos fc 2 — r cos fci — rcosfc 2 ). (A7) 


2. Electron-phonon coupling 

For this part, it is more convenient to write everything 
in terms of k\ and fc 2 explicitly. 

{—,k|Tj*f)|+,k) = t sin ki<7^ + — r cos fcicr^ -1- , 

(—, k||+, k) = t(sinfc 1 cr^ + + sin fc 2 crf + ), 


K 


hence, 

iv k = (-,k|r 1 ( f ) |+,k)(+,k|r 1 ( 2 p) |- I k) 

t sin 2 fci (—*/i 3 i?k — (A k — h\) sin 0 k cos 0 k ) 

+ f sin fcr sin fc 2 (E%. sin 2 0 k + fc 2 cos 2 0 k ) 

— r cos k± sin ki(iE^hi — h^hz) 

— r cos fci sinfc 2 (—i£ , k /i 2 — h^h\) 


(A8) 


Thus the adiabatic curvature is 

8<»>(k) = -2I„,{^L 5 } 


t 

zeE 


(2 E k ) 

t sin 2 k\ (to — r cos k\ — r cos fc 2 ) 

+ rt cos fci sin 2 k\ — rt cos k\ sin 2 fc 2 

fci) 


t 2 (sin 2 fci (m — r cos fc 2 ) — r sin 2 fc 2 cos 
2^ 


and it can be made isotropic 

r(pAi ^ _ i2sin2 fc i( m - 2r cos fc 2 ) 


tv 

f 2 sin 2 fc 2 (to — 2r cos fci) 

4^3 ■ ( A9 ) 
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